An Efficient Time-Domain Method to Model 
Extreme-Mass-Ratio Inspirals 



Priscilla Canizares 1 and Carlos F. Sopuerta 2 

Institut de Ciencies de l'Espai (CSIC-IEEC), 

Facultat de Ciencies, Campus UAB, Torre C5 parells, Bellaterra, 08193 Barcelona, Spain 
E-mail: 1 pcm@ieec .uab. es , 2 sopuerta@ieec .uab . es 

Abstract. The gravitational-wave signals emitted by Extreme-Mass-Ratio Inspirals will be 
hidden in the instrumental LISA noise and the foreground noise produced by galactic binaries 
in the LISA band. Then, we need accurate gravitational-wave templates to extract these signals 
from the noise and obtain the relevant physical parameters. This means that in the modeling 
of these systems we have to take into account how the orbit of the stellar-mass compact object 
is modified by the action of its own gravitational field. This effect can be described as the 
action of a local force, the self-force. We present a time-domain technique to compute the self- 
force for geodesic eccentric orbits around a non-rotating massive black hole. To illustrate the 
method we have applied it to a testbed model consisting of scalar charged particle orbiting a 
non-dynamical black hole. A key feature of our method is that it does not introduce a small 
scale associated with the stellar-mass compact object. This is achieved by using a multidomain 
framework where the particle is located at the interface between two subdomains. In this way, 
we just have to evolve homogeneous wave-like equations with smooth solutions that have to 
be communicated across the subdomain boundaries using appropriate junction conditions. The 
numerical technique that we use to implement this scheme is the pseudospectral collocation 
method. We show the suitability of this technique for the modeling of Extreme-Mass-Ratio 
Inspirals and show that it can provide accurate results for the self- force. 



1. Motivation 

Thanks to the development of Gravitational Wave Astronomy, the possibility of studying and 
understanding new features of the Universe is becoming a reality. In this regard, there is an 
ongoing effort to learn about all the physical aspects of one of the main sources of Gravitational 
Waves (GW) for the future space-based Laser Interferometer Space Antenna (LISA)p]: Extreme- 
Mass- Ratio Inspirals (EMRIs). EMRIs are formed when a massive black hole (MBH) located at 
a galactic centre, with masses in the range M, = 10 4 — 1O 7 M , captures a stellar-mass compact 
object (SCO), with masses in the range m* = 1 — W 2 Mq (white dwarf, neutron star or a stellar 
black hole). The SCO inspirals towards the MBH following a long series of highly eccentric orbits 
that shrink due to the loss of energy and angular momentum through the emission of GWs. This 
inspiral is produced by backreaction effects related to the action of the SCO's gravitational field 
onto its own trajectory. In order to model EMRI systems and their GW emission in a way that we 
can produce useful GW templates for analyzing the data stream produced by LISA observations 
we need to consider the details of the gravitational backreaction responsible of the inspiral. This 
is a very challenging task for theorists. Nevertheless, the problem can be simplified due to the 
extreme mass ratios of these systems, in the range fj, = m*/M, ~ 10~ 7 — 10~ 3 . This allows us to 



describe them in the framework of BH perturbation where the SCO is modeled as an accelerated 
point-like mass in the MBH background and the gravitational backreaction is pictured as the 
action of a local force, the self-force. For the purposes of our work, we can further simplify the 
problem by studying an analogous EMRI system which consist on a scalar charged point particle, 
q, orbiting around a Schwarzschild MBH (see, e.g. |2J) that fixes the spacetime geometry. In this 
simplified model, the inspiral proceeds due to the emission of a scalar field, generated by the 
SCO motion and which affects the SCO trajectory through the action of a self-force given by the 
gradient of the scalar field: 

F» = q(g^+uV)(V^)\,, «" = ^T> (!) 

where 7 and denote the SCO's trajectory and unit four-velocity respectively. We use this 
simplified model as a test bed to illustrate the techniques that we have developed for self-force 
computations. 

In the case of non-rotating MBH, the spherical symmetry of the BH spacetime provides 
additional simplifications of the problem. In particular, the scalar field can be decomposed into 
scalar spherical harmonics, Q im (t,r), which are decoupled between them. On the other hand, 
the point-like description of the SCO produces divergences in the retarded scalar field and hence 
we need to regularize it. In order to do so, we use the mode sum regularization scheme |3], 
which provides an analytical expression for the singular contribution of the retarded field, <I> S , 
at the particle location. In this way, by subtracting it from the full retarded field, <3?, we obtain 
a smooth and differentiable field at the particle location, $ R = $ — <I> S . Then, the meaningful 
expression for the self-force is: 

F" = g(r + M V)(V/)| 7 . (2) 

Thus, what we need is a (numerical) technique to compute the full retarded field, <I>, and to use 
the mode sum scheme to obtain the regularized self-force. In what follows we report on work 
we have carried out recently on the development of a new accurate and efficient technique for 
self-force computations of eccentric EMRIs [I]. This work is the extension of a previous one [SI [6], 
where these techniques were introduced for the case of circular EMRIs. 

2. Modeling EMRIs using a multidomain framework: The 
Particle-without-Particle Formulation 

Our computational scheme consists in diving the computational domain in a number of 
subdomains in such a way that the SCO (the particle) is always located at the interface between 
two subdomains. This has two main advantages: (i) We avoid introducing a spatial scale to 
resolve the point particle, (ii) The equations for the scalar field, which nominally have singular 
source terms, become homogeneous equations at each subdomain. These equations, assuming 
that appropriate initial date is precribed, lead to smooth solutions. This fact translates into 
good convergence properties of the numerical method used to implement this method. 

Because of these properties we call this formalism the Particle-without-Particle (PwP) 
formulation. We evolve the individual wave-type equations at each subdomain by using 
time domain methods, which perform well for eccentric orbits. This technique was already 
implemented for the case of circular orbits and has also been presented at the previous LISA 
Symposium (5J |6] . We have recently extended the method to make computations also in the case 
of eccentric orbits [4J. In what follows we summarize the main results of this work. 

We start by describing the multidomain structure of our PwP formulation (see Figure[T]). Once 
we have expanded the scalar field in spherical harmonics, each mode satisfies an independent 
(not coupled to the other modes) 1 + 1 wave type equation of the Regge- Wheeler type (with the 



*2 h 



tl h 



<o I- 









1— - - H 




h- - — -1 




1— -1 

1— - H 




N 
tip 



I 



f 3 h 



*2 h 
ti h 



t 



1 


* 


- H 


I-- - 




H 


1 — 


t - 


H 



I I 



vS(*0 

4F — i - 1- -— i 



W 



Figure 1. Structure of the one-dimensional spatial domain for the circular case (left) and for the 
eccentric case (right). While in a circular orbit the particle's radial coordinate remains constant 
and the grid is static, in a generic orbit case the particle oscillates between pericenter (r peri ) 
and apocenter (r apo ) and the grid has to readapt continuously to maintain the particle at the 
interface between two subdomains. 



potential associated to a scalar field). Then, the spatial domain is one-dimensional: 0, = [r* , r*], 
where r* is the radial tortoise coordinate, rjj corresponds to the truncation in the direction to 
the MBH horizon and r* corresponds to the truncation in the direction to spatial infinity. In 
order to describe the subdomain communication let us consider for the moment a splitting of 
the computational domain into two regions or subdomains, one to the left and one to the right 
of the particle: r* G [rjj, r*] = [rjj, r*] U [r*, r*]. Each of these regions can be in turn divided into 

more subdomains (see Figure lj: U = (jf = i wliere = [ r *a,L, r *a, R ] wi th < R = <+i, L - In 
the eccentric case (3], we need to use a coordinate system in which the particle is comoving with 
the interface of two such domains, or in other words, the coordinates L and r* R associated 
with the particle must be time dependent. 

In this setup, we use reduction the wave-like equation to a first-order system by promoting 
the time and radial derivatives to new variables. Then, these variables, U, can also be split as 
follows: 

U(t,r*) = U_(t,r*)e(r* p (t)-r*)+U + (t,r*)e(r*-r;(t)), (3) 

where G denotes the Heaviside step function. We can define the jump across the particle's 
location in an arbitrary quantity as: [A] = lim A. (t, r*) — lim \_(t,r*). We can obtain 

expressions for the jumps by inserting Q into the field equations (see |7] and jl]). Using 
these junction conditions we communicate the solutions of the homogenous equations at each 
subdomain. 

The numerical method that we use to solve our equations is the PseudoSpectral Collocation 
(PSC) method |8]. We discretize each subdomain by using a Lobatto-Chebyshev grid. Then, 
in the PSC method the variables have a spectral representation in terms of an expansion in 
Chebyshev polynomials, {T n (X)} [X G [—1, 1], n = 0, . . . ,N), and a physical representation in 
terms of the values of the variables at the collocation points: 

N N 

U N (t,r*) = Y, a n(t)T n (X(t,r*)) = £ U^t) C^X^t, r*)) , (4) 
n=0 i=0 

where X(t,r*) is the mapping between the spectral and physical domains (the time dependence 
appears only in the eccentric case and is the way in which we keep the particle at the 



interface between two subdomains), a n are the spectral coefficients, and C i are cardinal functions 
associated to the Lobatto-Chebyshev grid [8]. One can change from one representation to the other 
by means either of matrices or fast Fourier transforms (as we do in our implementation). An 
important feature of the PSC method is that it provides exponential convergence for smooth 
functions, which is the case of our solutions after applying the PwP formulation. This is 
illustrated in Figure |2j where we show some convergence plots. 
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Figure 2. Truncation error estimated from the absolute value of the last spectral coefficient, a N , 
of the variable ijr m = r<& im . We show the results for the harmonic modes (£,m) = (2,2) (left 
column) and (£,m) = (20,0) (right column) and for the orbital parameters (e,p) = (0.1,7.0) 
(top row) and (e,p) = (0.5,8.0) (bottom row), where here e denotes the eccentricity and p 
the semilatus rectum. We observe spectral convergence (straight line) until we reach machine 
roundoff error (plateau). The data for these plots has been obtained from the subdomain at the 
right of the particle. 

Once we have performed the spatial discretization using the PSC method, we obtain a set of 
ordinary differential equations that can be evolved using the Method of Lines. The particular 
implementation that we use involves a Runge-Kutta 4 solver. The evolution must include 
the communication between subdomains that we have discussed above. We use two different 
numerical techniques to communicate the subdomains: (i) The penalty method, where the 
system is driven dynamically to satify the junction conditions, and (ii) the direct communication 
of characteristic fields, where the junction conditions are imposed by communicating the 
characteristic fields of the first-order system of partial differential equations. To illustrate the 
ability of our method to resolve the field and its derivatives around the particle we show some 
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Figure 3. Snapshots of the evolution of the variables ip im (left), <j/ m = dtip irn (center), and 
ip im = d r tip im (right), for the mode I = m = 2. The orbital parameters are (e,p) = (0.5, 8). 

snapshots of the evolution in Figure [3] 
3. Some Results and Discussion 

Until now we have introduced the foundations of our method and the particular numerical 
techniques that we use in order to implement it. We have shown that the PwP formulation lead 
to smooth solutions at each computational subdomain, which in turn exploits the exponential 
convergence of the PSC method as shown in Figure [2] We have also shown how well this method 
can resolve the field and its derivatives around the particle location (see Figure [3]). In particular, 
we can see how the jumps in the time and radial derivatives of the scalar field modes & m are 
well resolved (see the inset plots in Figure [3] where we zoom in the area near the particle) . 

Finally, we have computed the values of the gradient of the regularized components of the 
gradient of the scalar field, V Q $ fl , which contain information equivalent to the self-force (see 
Eq. ([2])), which is the projection orthogonal to the particle four-velocity of V a & R . In Figure [4] 
we show the evolution of the components of the regularized field obtained by truncating the 
sum over harmonic modes in such a way that we only include modes with £ ^ £max — 

20. This 

figure shows the evolution of the regularized field for two different orbits, which in terms of 
the eccentricity and semilatus rectum are: (i) (e,p) = (0.1,7.0) and (ii) (e,p) = (0.5,8.0). In 
order to quote some numerical results, we have computed the values of the self- force components 
using the method of the comunication of the characteristic field described in Sec. [2] We obtain: 
$f = 2.2824 x 10~ 4 (g/M. 2 ), $^ = 8.5836 x 10~ 5 (q/M?) and = -3.7129 x 10~ 3 (q/M,) 
for the eccentric orbit (i) and $f = 6.4751 x 10~ 6 (q/M?), = -5.3519 x 10~ 6 (q/M?) and 
= 2.8778 x 10~ 4 (q/M,) for the eccentric orbit (ii). 

All our calculations have used between 8 and 39 subdomains and 50 collocation points 
per domain. The average time for a full self-force calculation (which involves the calculation 
of 231 harmonic modes) in a computer with two Quad-Core Intel Xeon processors at 2.8 
GHz is always in the range 20-30 minutes. These calculations can be further optimized by 
distributing the subdomains and collocation points so that the resolution is adapted further 
to our physical problem, and is this is the subject of ongoing work. The calculations can be 
easily parallelized, either by spreading the computational calculations for individual modes or 
for individual subdomains. In the last case the parallelization is not trivial as we need to pass 
the relevant information for subdomain communication. 

Looking to the future, we are currently working on extending these techniques to the 
gravitational case, where the challenge comes from the fact that each harmonic mode is described 
by a set of coupled 1 + 1 wave type equations, but where the PwP formulation can be used as it 
has been described here. 
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Figure 4. Evolution of the components of the gradient of the regularized field, V a $ R , for a 
scalar charged particle in eccentric orbits around a non-rotating MBH. From left to right, the 
orbital parameters of the orbits are: (i) (e,p) = (0.1,7.0) and (ii) (e,p) = (0.5,8.0). For each 

orbit (frame), the solid line represents the evolution of the dimensionless time component, * <I>f ; 

M 2 

the dashed line represents the evolution of the dimensionless radial component, ■ — - and the 
dot-dashed line represents the evolution of the dimensionless azimuthal component, 
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